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We study the magnetization of square and hexagonal graphene dots. It is shown that two classes of hexagonal 
dots have a second order phase transition at a critical Hubbard energy U, whose value is similar to the one in 
bulk graphene, albeit the dots do not have a density of states proportional to the absolute value of the energy, 
relatively to the Dirac point. Furthermore, we show that a particular class of hexagonal dots having zig-zag 
edges, does not exhibit zero energy edge states. We also study the effect of uniaxial strain on the evolution of 
the magnetization of square dots, and find that the overall effect is an enhancement of magnetization with strain. 
The enhancement can be as large as 100% for strain of the order of 20%. Additionally, stress induces a spatial 
displacement of the magnetization over the dot, moving it from the zig-zag to the armchair edges. 

PACS numbers: 73.20.-r, 73.20.At, 73.21.Ac, 73.22.-f, 73.22.Gk, 81.05.Uw 



I. INTRODUCTION 

Nowadays, the terms wonder material and graphene 
dreams^ ^ freq uently accompany the description of the un- 
usual electronic!^ thermal^, and mechanical properties of 
graphene^ ^ One of the most promising graphene dreams 
is its application to a new generation of nanoelectronic 
devices^ ^. To that effect, a number of systems have already 
been experimentally investigated, namely: single-electron 
tra nsistor ^, quantum interference devices and graphene 
(Jq^Misj jjjg presence of Coulomb oscillations in graphe ne 
quantum dots was also identified by different group s^^HEII 

Theoretically, the first investigations in this context focused 
on the transport properties of short (and wide) ribbons 
For long graphene ribbons it was shown that the low bias 
current flowing through the bulk of the ribbon is very ro- 
bust with respect to a variety of constriction geometries and 
edge defects, a result also confirmed for disordered armchair 
nanoribbonJI^. As in the case of short ribbons evanescent 
waves were seen to play an important role in the electronic 
transport through graphene quantum dots^^. The role of mag- 
netic fields in the electronic levels of graphene quantum dots 
has been investigated by several authors. Of particular interest 
for transport properties is the fact that optical properties can 
be tuned by the size and edge type of the dot^^. The shape 
and symmetry^^^ of the dots also play an important role on 
energy level statistics and charge density. For the special case 
of triangular quantum dots^^, the existence of "ghost states" 
was revealed, when these dots have armchair edges, whereas 
for triangular dots with zigzag edges the well known surface 
states are present. Of particular importance was the demon- 
stration of large insensitivity of the electronic structure to the 
edge roughness^^. 

The main motivation for research in graphene quantum dots 
and ribbons is related to the need of producing a graphene- 
based system with an energy gap, which is not present in 
bulk graphene. This fact is a recognized shortcoming of 
bulk graphene, in what concerns applications relying on cur- 
rent electronic operation. Gaps can be induced by electronic 
quantum confinement in narrow armchair ribbons^^, a result 
confirmed by ab-initio calculations^^, and experimentall^^"^. 



First principle calculations further show that zig-zag ribbons 
can support magnetic ground states which leads to a gapped 
spectrunP3. Spin polarized ground states are also present in 
small graphene derivative molecular systems'^. This find- 
ing opens a new line of research: the study of spin polarized 
ground states of graphene quantum dots of different geome- 
tries. In both single^^ and bilayer^^ zig-zag ribbons, it was 
found that opposite edges align antiferromagnetically, with 
the magnetization rapidly decaying towards the bulk of the 
ribbon. Hartree-Fock (HF) calculations are specially suited 
for this type of study, since one can study the effect of dif- 
ferent values of the Coulomb interaction on the magnetic 
structure. An interesting effect emerges from such studieJ^ 
graphene ribbons have no critical Hubbard interaction, U, and 
thus the HF ground state is always magnetic. This result is at 
odds wit h the b ehavior of bulk graphene^^ 29 so 3 1 bilayer 
This richness of different behaviors suggests 
studying the formation of magnetic ground states in graphene 
dots, a line a research we carry on in this paper from an HF 
point of view. 

Interest in magnetism of sp^ carbon systems was greatly 
spurred by experiments with proton-irradiated graphite^^, and 
with the experimental evidence that the measured magnetism 
might stem from 7r-orbital physics alone^^. Proton irradia- 
tion induces spatially disordered vacancies in the system'^^EU 
The magnetism found experimentally is supported theoreti- 
cally by Hartree-Fock and ab-initio studies"^^. Recent exper- 
imental developments addressed the intrinsic ferromagnetism 
in HOPG graphite, originating from the naturally occurring 
grain boundaries, where zig-zag edges develop and local mag- 
netic moments are formed. Typical hysteretical curves of a 
ferromagnetic material are seen in a temperature range from 
5 K up to 300 K^^ 

Hence, disorder, such as those line defects studied recently, 
is a possible route for ferromagnetism in carbon-based mate- 
rials. However, disorder is not a necessary condition for mag- 
netism in sp^ systems. It is by now well established theoret- 
ically that graphene systems with zig-zag edges can support 
magnetic moments and, in a system with perfect edges, this 
leads to magnetic ground states. An argument widely used 
against this result is based on the fact that the atoms at the 



2 



edges are, essentially, gigantic free radicals, which would be 
impossible to realize in a true graphene system. This argument 
ignores, however, the fact that such gigantic free radicals can 
be chemically passivated with other chemical species, notably 
hydrogen. It is found ab-initio that the long range mag netic 



order is robust, even under passivation of the edges^, con- 
firming early predictions^. 

In small graphene structures (triangles and hexagons), mag- 
netism has been thoroughly investigated^^ both using ab-initio 
and Hartree-Fock methods, but generally for small dot sizes. 
The interplay between transport and magnetism has also been 
addressed"^^, as well as magnetism induced by vacancies'^ ^ 
(which can be seen as a three site zig-zag edge). 

Another topic of experimental research that has recently 
seen a considerable upsurge, is the study of interplay between 
the mechanical properties of graphene and its electronic struc- 
ture. The motivation for these studies is the possibility of tai- 
loring the transport properties of graphene by means of ex- 
ternally induced strain^. Naturally related is the question of 
how can the above mentioned magnetic properties of graphene 
ribbons and dots be modified by external stress. In a previous 
work"^^, some of the present authors showed that the electronic 
spectrum of graphene can be strongly modified by external 
stress. In particular, stress along the zig-zag edges of the sys- 
tem might eventually lead to the opening of a gap at large 
deformations. In addition to studying the magnetic properties 
of quantum dots in equilibrium, here we will also address how 
their magnetization is affected by external stress. 

Our main findings can be summarized as follows. The ex- 
istence of magnetic ground states in graphene dots of nano to 
mesoscopic sizes depends on their geometry, and not only on 
the existence of zig-zag portions along their edges. Within a 
Hartree-Fock framework, the existence (or not) of a minimum 
on-site Coulomb repulsion, Uc, for the onset of magnetism 
depends critically on the dot geometry and symmetry. When 
strain is applied, the nearest-neighbor hopping integrals are 
naturally modified. This leads to a modification of the local 
magnetic moments found in the ground state, in a way which 
is much stronger than one would expect just by calculating the 
isotropic renormalization of the critical Coulomb repulsion U. 
Our results show that magnetism is enhanced under uniaxial 
strain, and causes a reduction of Uc, for the dots which exibith 
finite Uc. Moreover, we find that, under strain, the local mag- 
netic moments associated with zig-zag edges in rectangular 
dots can drift from the zig-zag to the armchair edges. 

The paper is organized as follows: in Sec. |Il|we introduce 
our theoretical model, and discuss the relevance of several 
Coulomb terms in defining an effective Coulomb interaction 
U. A discussion of the appropriate value of U for graphene 
ensues. In Sec. |III| we study the magnetization of different 



types of square and hexagonal graphene dots. In Sec. [lV|the 
role of strain on the magnetization of graphene dots is consid- 
ered. Our main results are discussed in Sec. (Vl 



II. MODEL 

Our study of magnetism in graphene quantum dots and an- 
tidots relies on the Hubbard model with on-site interaction, 
an approach used by other authors in the study of graphene 
ribbons^^ . Since dots have no translation symmetry, the prob- 
lem is solved in real space. To that end, we need to set up the 
Hamiltonian in a matrix form, which requires a convenient al- 
gorithm to build such a matrix for the different type of dots. In 
what follows we describe the model, together with the physi- 
cally relevant values of the on-site Coulomb interaction. 

The study of magnetism in condensed matter physics is 
traditionally, and frequently, based upon the Hubbard model, 
which can be written as 



H = Hq -\- Hu 



(1) 
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4(r)6<,(r + ,5) + H.c., (2) 

Hu = U^a\{r)ai;{r)a\{r)ai{r) , 

r 

+ uY,b\{r)b^{r)b\{r)bi{r). (3) 



For graphene the hopping integral is t :^ 2.7 eV (used as the 
energy unit in this work), U is the on-site Coulomb repulsion 
energy, and aj.(r) [b]^{r)] is the electronic creation operator 
at site A [site B] of the unit cell r in the honeycomb lattice 
(see Fig. [T]). 

The Coulomb term is treated at the mean-field leveP^l^by 
making the replacement of the quartic interaction by 
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U ^al{r)a^{r){o}_^{r)a_^{r)) , 

r,cr 

uY,bl{r)b,{r){bl,{r)b.,{r)) , (4) 



such that, when cr =T,i, we have —cr =i, j. After this 
transformation, the quantum problem becomes bi-linear in the 
electronic operators, and can be solved by diagonalization of 
two matrices of dimension D x D, where D is the total num- 
ber of lattice sites in the dot. The electronic density has to 
be determined self-consistently and the mean-field equations 
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na,a{r) = {al{r)a^{r)) , 
n,,,{r) = {bt{r)b,{r)) . 



(5) 
(6) 



Here na,cr(^) and nb^cr{r) are the mean electronic densities of 
spin cr at the A and B sites of the unit cell r, respectively. The 
wave function of the system corresponding to an energy Ex^a 
is labeled by the quantum number A, having the explicit form 

\^x,a) = ^ AA,a(r)|a,r) + BxAr)\b,r) , (7) 



where |a,r) and |6, r) are lattice-position basis-states. The 
mean field equations ^ and ^ are determined as function 
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FIG. 1 : (Color online) Illustration of the honeycomb lattice with the 
A and B sublattices, the lattice vectors (5^ (i = 1, 2, 3), the primitive 
vectors a and 6, and the hoppings ti (i = 1, 2, 3) used in Sec. |IV] 
The abscissas are along the zig-zag edge (horizontally in the figure), 
and ao is the equilibrium carbon-carbon distance. 



of the AA,cr(^) and Bx^^{r) coefficients, defined in the wave 
function (|7]), according to 



na,a(r) = ^|AA,.(r)p/(^A,a), 

A 

n5,.(r) = ^|5A,a(r)|'/(^A,a), 



(8) 
(9) 



where f{x) = (1 + g/3(^E-/i) j-i^ ^ ^j^^ chemical potential 
and (3 = l/(/c^T), with T the temperature. The problem 
has to be solved numerically. We start with a trial solution 
for na^o-(^) and n5^cr(^); then the Hamiltonian is diagonalized 
and new values for na,cr(^) and n5,cr(^) are computed; the 
procedure is iterated a number of times until convergence is 
reached. 

As mentioned, the conventional treatment of magnetism in 
graphite and graphene at the Hartree-Fock level includes only 
the effect of the on-site Coulomb interaction U. We now dis- 
cuss the importance of more general interactions^^. We first 
note that, at the mean-field level, a nearest neighbor Coulomb 
interaction does not contribute to the existence of a ferromag- 
netic phase in the case of a system with translational invari- 
ance and a single orbital per unit cell^^. If in graphene we 
consider a Coulomb term of the form 

Hy=V al{r)a„{r)bl{r + 5)b„,{r + 5), (10) 

it remains true that such interaction will not contribute (at the 
Hartree-Fock level) to the existence of a magnetic ground state 
in the thermodynamic limit. The situation is different, though, 
in a system without translational invariance, since the spin 
density in neighboring carbon atoms is not necessarily equal. 
This is of special relevance near the edges of the system. 
The mean field Hamiltonian has the form 

Hy ^ H^^ = vY,al{r)a,{r)n^{r) + (a ^ 6) , (11) 



h operators interchanged, and 

5,ct' 

is the average density at the B neighbor carbon atoms of a 
given A atom at position r. The terms Hu and Hy are direct 
Coulomb interactions. An exchange term can also be included 
in the Hamiltonian, having the form 



Hj = -Y,S,{r)-S,{r + 5)., 

r,d 



(13) 



with Sa{r) [Si){r)] the electronic spin operator of an electron 
at site r of the sub-lattice A [B] . In this case, the mean field 
Hamiltonian is 



Hj^Hy^ = -J2 4(r)a.(r)S,(r) + (a ^ b) , (14) 

r,a 

with 



where (a ^ b) in Eq. ( 1 1 ) is a shorthand notation for a term 



with the same form as the first, but with the role of the a and 



= ^a^iblir + S)b^^{r + S)) , (15) 



where S6(r) is the average spin density at the B neighbor of 
a given A atom at position r, and a takes the values ±1 when 
used as a multiplicative factor. 

We shall assume that the leading overall effect of these three 
interactions can be captured by a renormalized Hubbard inter- 
action, U, in the mean field calculations. Thus the value of 
U should reflect this effective interaction, rather than the bare 
on-site Coulomb repulsion in graphene. 

We now proceed to study the ground state of dots and their 
magnetization as a function of U. A natural question imme- 
diately arises: what value of effective U should one take to 
be consistent with the magnitude of the real Coulomb interac- 
tions in the material. For the benzene molecule U was seen 
to be as large as 16 eV"^^. In a recent study of magnetism in 
disordered graphene and irradiated graphite"^^, the value of U 
was considered to be in the interval 3-3.5 eV, based on the 
value accepted for ^ra^^-polyacetylene, a one-dimensional bi- 
partite sp^ carbon system (although this value of U for ^ra^^- 
polyacetylene has been subject to controversy^^). Other two 
recent studieJ^^^ took U =2 eV and U = 3.85 eV, values 
that reproduce the LDA gap in graphene ribbon^ and the 
HOMO-LUMO gap in small graphene based structures. We 
shall consider below /7 = 2 eV and U = 3.5 eV as reference 
values in our calculations. 



III. A (iV) AND MAGNETIZATION 

The results for the magnetization of graphene dots depend 
on the type of edges present. Generally speaking, one expects 
to see larger magnetization close to zig-zag edges, where the 
existence of localized states satisfies a spatial Stoner crite- 
rion for finite values of U^^. The existence of such type of 
states is shown to be related to lattices with an odd number of 
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sites^^ . For models with sub-lattice symmetry, as is the case 
of graphene, the number of zero energy modes is determined 
from the difference |7V^ — Nb\, were Na and Nb is the num- 
ber of sites in sub-lattice A and B respectivel>^^^^. 

In our calculations we use relatively small dots. This choice 
is justified because there are almost no visible finite- size ef- 
fects, as we explicitly show below by studying dots of differ- 
ent sizes. Additionally, our choice is also justified from an 
experimental point of view, since it recently became possible 
to cleave graphene crystalites down to one-dimensional chains 
by irradiation inside a transmission electron microscope [see 
Fig. 2 of Ref. ^|. With such new experimental methods, 
tailoring dots of any possible size and shape seems now quite 
within reach. 

It is useful for latter use to introduce the quantity A(A^), as 
the energy interval between the highest hole state and the first 
particle one, for the system without interactions (U = 0): 

A(iV) = i?C:if -i?httst, (16) 

where N is the total number of atoms in the dot. We consider 
two types of dots, with square and hexagonal shapes, and also 
the case of a dot with two non-connected regions (some times 
referred to as an anti-dot). 

We start with the study of hexagonal dots. There are hexag- 
onal dots with different symmetries and different types of 
edges: 

1. dots with Dq symmetry, having only armchair edges 
(see Fig.[2]a) ). 

2. dots with Dq symmetry, having armchair and zig-zag 
edges (see Fig. [3]b) ). 

3. dots with Ds symmetry, having armchair and zig-zag 
edges (see Fig|2]c) ). 

The first type of dot defined above shows that it is possible 
to have dots without zig-zag edges, no matter how large they 
are, and therefore the physics associated with zig-zag edges 
should not be present. This type of Dq dot, when very large, is 
almost equivalent to the bulk system, having the full symme- 
try of the honeycomb lattice and therefore showing a second 
order phase phase transition at a (mean-field) critical Hubbard 
interaction, Uc, given by 

Uc^2.23t, (17) 

as shown in Fig. |4](HEX2-type). The same holds true for the 
Dq hex 1 -type of hexagons, but with a smaller value of Uc 
(smaller than 2). The dependence of the maximum value of 
magnetization as a function of U for the two Dq hexagons is 
plotted in Fig. |4] There we see that the critical U is close to 
that given by Eq. ([17]), without any noticeable variation with 
the size, L, of the hexagon. In Fig. |4]the reference values for 
U discussed at the end of Sec. [lljare represented as vertical 
dashed lines. Clearly, the magnetic transition is well above 
those reference values for U, meaning that this type of dots, if 
experimentally fabricated, should exhibit no magnetic order. 




FIG. 2: (Color online) Types of hexagonal and square dots studied 
in this work: (a) hexagon with Dq symmetry and no zig-zag edges 
(termed HEX2 in the figures below); (b) hexagon with Dq symmetry 
and zig-zag edges (termed HEXl in the figures below); (c) hexagon 
with Ds symmetry with zig-zag edges; (d) an anti-dot with the ex- 
ternal and internal boundaries made of Dq symmetry with no zig- 
zag edges; (e) square where the vertices are of zig-zag type (termed 
SQRl); (f) square where the vertices are of armchair type (termed 
SQR2). The size of the figures is characterized by the number L of 
carbon atom horizontal-lines (zig-zag type of lines); for example, in 
panel (a) one has L = 14, and in panels (e) and (f) one has in both 
cases L = 10. 

We note that the two Dq hexagons have finite A (A/") values, 
which vary as a power law with N, as shown in Fig. |5] For the 
HEXl- and EXH2-types of hexagons we numerically extract: 

A(7V) 1.717V-^•^^ (18) 
A(7V) 1.7bN-^-^^ , (19) 

respectively. The exponent in the above power laws is essen- 
tially equal to ^, and, therefore, reflects the finite-size quan- 
tization of the electronic spectrum. For square dots, on the 
other hand, we find that A (A/") vanishes much rapidly as A^ 
increases, reflecting the formation of edge states at nearly zero 
energy: for small systems, the edge states from opposite sides 
of the square dot hybridize, and the otherwise zero energy 
states for the semi-infinite system split in energy. As the width 
of the dot increases the hybridization is strongly suppressed 
and zero energy levels develop. 

The finiteness of A (A") for the hexagons correlates with 
the finite value of Uc seen in Fig. |4] On the other hand, the 
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value Uc — 2.23t previously obtained in the literature^ was 
determined using the fact that the density of states of bulk 
graphene is proportional to the absolute value of the energy 
relatively to the Dirac point, being zero for a half-filled sys- 
tem. These two results - for Dq hexagonal clusters and the 



bulk system - means that the value of Uc ( ^7] ) is not exclu- 
sively determined by the vanishing nature of the density of 
states at the Dirac point of bulk graphene. On the other hand, 
the two Dq hexagons show different values of Uc, which can 
only be interpreted as a boundary effect, determined by the 
different nature of their edges. It is worth noticing that the 
hexagons of type HEXl, having zig-zag terminations (defin- 
ing a figure with Dq symmetry) do not develop zero energy 
states, leading, therefore, to the finiteness of Uc. The hexago- 
nal dot of Ds symmetry shows a behavior for A{N) identical 
to that found for the squares and, as a consequence, there is no 
finite value of Uc'. the system is magnetic for any arbitrarily 
small value of U. 

The mean-field values of Uc determined for the Dq 
hexagons will be modified by quantum fluctuations. The ef- 
fect of quantum fluctuations amounts in general to shifting the 
Hartree-Fock Uc to higher values^^^^^^. In the case of small 
graphene-based nanodots, such as bisanthrenJ^(C28Hi4), the 
Hubbard Coulomb interaction may be larger than that as- 
sumed for macroscopic sp^ carbon systems. This hypothesis 
is based on the value U ~ 16 eV computed for benzene'^^ 
Given this value and our current results, there is a real pos- 
sibility of having magnetic ground states in small hexagonal 
systems with Dq symmetry. 

In what concerns the relation between magnetism and edge 
structure, we see that dots with Ds and square symmetry have 
zig-zag edges, and this leads to finite magnetization for any fi- 
nite U. Magnetization is maximal at, and close to, the zig-zag 
edges, and fades rapidly as one progresses towards the bulk 
of the dot. It is worth noticing that, for HEXl -type hexagonal 
dots, there are six external zig-zag boundaries, but its spec- 
trum does not present zero-energy eigen values. In Fig. |3] 
we present particular cases of the spatial distribution of mag- 
netization in the different dots considered in Fig. [2] For the 
hexagons of type HEX2 the magnetization is homogeneous 
over the boundary and, as soon as /7 > /7c, it develops from 
the boundary of the hexagon toward the bulk. For HEXl 
hexagons the maximum of magnetization develops at the zig- 
zag vertices, but again only for U > Uc. In the case of the 
hexagon with Ds symmetry the development of the magneti- 
zation follows the pattern of that found for HEXl dots, but it 
is finite for any finite U value. The antidot case (panel (d) of 
Fig. |3]) can be considered a simple case of a disordered sys- 
tem, since all symmetries are broken. In this case the magneti- 
zation develops preferentially at the internal edges connected 
to the bulk of the system, with formation of a shadow region 
at bottom left of the anti dot where no magnetization is seen 
(for that particular value of U). This behavior can be under- 
stood since the internal boundary plays, in the anti-dot, the 
same role as the external boundary of the equivalent dot. It 
is worth noticing that the anti-dot, being made of HEXl and 
HEX2 hexagons has a critical Hubbard interaction, which is 
controlled by the Uc value of the HEXl hexagon. 



a) 































FIG. 3: (Color online) Illustration of the local magnetization for dots 
of the same type illustrated in Fig. [2] but with a much larger num- 
ber of atoms, thus avoiding finite-size effects. Upright triangles refer 
to positive magnetization and the down ones to negative values. The 
open triangles refer to the points where the magnetization has it max- 
imum value. In panel (a) we have a HEX2 hexagon, with U = 2.3; 
for (b) and (c) hexagons [/ = 2. In panel (d) we show an anti-dot, 
where the external boundary is from a EHX2 hexagon and the inter- 
nal one is from a HEXl hexagon. Panels (e) and (f) are of type SQRl 
and SQR2 respectively. 



In Fig. |4] we depict the dependence of the maximum of 
magnetization (mmax) with U. We see that for the square dots 
of both types considered in Fig. [2] the magnetization is finite 
down to arbitrarily small values of U. This can be correlated 
with the correspondingly small values of A(A^), shown if Fig. 
|5j and the behavior of the DOS at the Fermi level. As to the 
hexagonal dots (see Fig. |4]), we also see that the existence of a 
finite A(A^) is associated with the existence of a finite Uc. In 
other words, some geometries have finite density of states at 
E = 0, p{0) 7^ 0, whereas others do not. In the case p(0) = 
and as A/" ^ oo the behavior of the system is essentially that 
of the bulk case up to finite-size corrections. On the other 
case, with p(0) 7^ for finite N, the magnetic behavior is 
different and there is finite magnetization for any value of U. 
We then understand the result obtained by Sorella and Tosatti 
C^) as the limiting case of A' ^ 00 with p(0) = for any 
finite A^. 

A relevant quantity to compute is the energy difference be- 
tween the paramagnetic and ferromagnetic ground states, de- 
fined as 



AE = E^ 



Efe: 



(20) 



where £^para and £^ferro are the ground state energies of the 
paramagnetic and ferromagnetic ground states. The value of 
AE is intimately related to the (mean-field) temperature at 
which such a magnetic ground state could be observed. In 
Fig. |6]the value of AE is given in Kelvin for both hexagons 
and squares. As discussed previously, these squares magne- 
tize for any finite U and therefore we could, in principle, ob- 
serve magnetism with the values of U expected for graphene 
(dashed lines in Fig. [6]). The mean-field critical temperature is 
relatively high, between 10 K and 30 K, in the interval for the 
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FIG. 4: (Color online) Variation of the magnetization as function of 
U, for squares and hexagons of different types. The vertical dashed 
lines refer to the values of U used in Refs. I43l51l (see text in Sec. |ll| 
for a discussion about these choices). 
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FIG. 5: (Color online) Variation of quantity A(N) with the num- 
ber of atoms for square and hexagonal dots. From this figure alone 
we can understand that the spectrum of hexagons of type HEXl and 
HEX2 behaves exactly in the same way, not showing the develop- 
ment of zero energy edge states, a consequence of the D6 symmetry 
alone. 



lower and higher expected U in graphene, but not comparable 
to the room temperature magnetism observed in graphitd^. 



IV. MAGNETISM AND STRAIN 

Strain in graphene in now an active topic of experimental 
research. It was shown that some amount of strain can be 
induced either by deposition of oxide capping layers^ or by 
mechanical methods^^. The amount of strain can be deter- 
mined by monitoring the blu#^or recP^ shifts of the G and 2D 
Raman peaks of graphene. This method is a straightforward 
extension of related studies used in graphite nano-fibers^2 
Strain has also obvious consequences on the electronic and 
heat transport, producing metal- semiconductor transitions, as 
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FIG. 6: (Color online) Variation of AE with the Coulomb interaction 
U/t, for hexagons of the type HEXl, and squares of the type SQRl. 



in carbon nanotubes'^, or transport anisotropy in graphenJ^. 
These effects are due to changes in the band-structure of the 
materials as a consequence of the modification of inter-atomic 
distances, which in turn implies a change of the electronic 
hopping parameters. To our best knowledge, the first correct 
studies of strain effects on the bandstructure of graphene were 
undertaken in Ref . 47 62 

For hexagonal systems the relation between stress, cr^j, and 



strain, Uij (i,j = x^y,z), read^ 



^yy 



"yz 

"^zx 
'^xy 



Su 

Sis 






5*12 

^11 
Sis 









Sis 













^xx 


Sis 













^yy 


Sss 




























^yz 








6*44 







C^zx 











2{Sii-Si2)_ 




_ ^xy _ 

(21) 



where the elements Sij (here i^j = 1,2,3,4) are termed 
compliance constants. For the case of graphene under uni- 
axial tensile strain, the relation between stress and strain is 



'XX — Slid XX 1 
^yy — Si20'xx i 



(22) 
(23) 



meaning that graphene behaves as an isotropic elastic 
medium. We shall consider two cases of stress: applied along 
the zig-zag edges {ZZ), and applied along the armchair {AC) 
edges. In these two cases, the absolute values of the next- 
nearest-neighbor vectors 5i change aJ^ 



l^i.sl = 1 

|<52| = 1 
for the ZZ case, and 

l^l,3| = 1 

l^2| = 1 



1 



4 



1 3 

-e - -ev , 
4 4 



(24) 
(25) 

(26) 
(27) 



for the AC case, where e = ^ncr is the amount of longitudi- 
nal strain, and v = —S12/S11 is the Poisson ratio. The two 
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FIG. 7: (Color online) Representation of the effect of stress on the 
length of the nearest neighbors carbon atoms. On the left we de- 
pict the ZZ case, and, on the right, the AC one. For the ZZ case, 
the stress a induces an increase of the hopping associated with the 
vertical bond, due to the Poisson effect. The hoppings associated 
with bonds 1 and 3 are reduced, and the system tends to dimerize 
at large deformations. For tension along AC, a is oriented along 
bond 2. In this case all hoppings are reduced, but the one associated 
with bond 2 decreases more than the other two, leading to a set of 
quasi one-dimensional chains. The quantitative change of the hop- 
pings upon stress was studied quantitatively using ab-initio methods 
inRef. [62] . 

cases correspond to two different physical situations, which 
can be understood in the case of extreme deformations using 
a simple picture: in the ZZ case the system tends to dimer- 
ize, since |^i,3| lengthen, and 1^21 shortens; in the AC case, 
all three distances lengthen, but 1^21 lengthens more, which 
can be construed as a tendency for the formation of quasi 
one-dimensional structures. A real space picture representing 
these two situations can be seen in Fig. |7] 

In Fig. [5] we present results regarding the effect of strain on 
the magnetization of the dots. For the purpose of illustration 
we consider dots of smaller size, but the results of Fig. [4]guar- 
antee that we should have negligible finite-size effects. The 
global effect of tensile stress on the edge magnetization of the 
dot is an increase of its magnitude, independently of whether 
we consider stress along ZZ or AC directions. In quanti- 
tative terms, the magnetization increase is more pronounced 
when stress is applied in the ZZ configuration than on the 
AC one. In experimental terms, the prediction is that mag- 
netism in graphene-based systems should be easier to detect 
when the material in under stress. 

Figure [9] shows the explicit variation of Uc with strain, for 
the hexagons of type HEX, which have a finite Uc- In the same 
figure, on the bottom row, the effect of e on the maximum of 
the magnetization, mmax, is also represented for squares of 
the type SQR2 (which have no critical U). We observe an 
increase of mmax as e increases. The dependence of mmax 
on e is not the same for the ZZ and AC cases at the same 
U value. Stress along the zig-zag edge is most effective at 
producing an enhancement of mmax- This behavior can be 
understood on the basis of the qualitative physical picture de- 
scribed in Fig. [tJ stress along zig-zag edges tends to pro- 



HEX2, L=8 SQR2, L=10 




U/t 



FIG. 8: (Color online) Variation of the magnetization as function 
of the Hubbard interaction [/, for different strain values e. The top 
panels refer to stress along the ZZ edge and the other to stress along 
the AC edge. The Poisson ratio used was that of a PET substrate 
(u ~ 0.3), considered in study of the Raman red shift of the G and 
2D peaks of graphend^. The vertical dashed lines refer to the values 
of U used in Refs. 14315 II (see text in Sec. |ll|for a discussion about 
these choices). 



ZZ AC 




FIG. 9: (Color online) Variation of the critical value, Uc, with the 
amount of strain, for the cases ZZ (top left) and AC (top right). In 
both cases we used HEX2-type hexagons with L = 8. In the bottom 
row we have the dependence of the maximum value of magnetiza- 
tion, mmax, on the amount of strain, using U/t = 2, and for the 
cases ZZ (down left) and AC (down right). In both cases we used 
SQR2-type squares, with L = 10. 



duce dimmers weakly coupled between them, which favors 
the magnetic state at those tightly bound atoms. 

For hexagons with a finite Fig. [5] shows that the over- 
all effect of strain along both the ZZ and AC cases is to re- 
duce the value of which for large £ obeys Uc <^ t. At 
first sight this result may seem easy to understand: the value 
of U cannot change with stress because it is a local (on-site) 
property^^. The hopping, on the other hand, depends strongly 
on the inter-atomic distance, and hence on the external stress. 
Since the result for Uc in the bulk system, Eq. ([17]), is bound to 



FIG. 10: (Color online) The top row shows the variation of the near- 
est neighbor hoppings ti,2,3 and the average hopping, (t) under uni- 
axial strain, with strain applied along the ZZ (left) and AC direc- 
tions (right). The bottom left panel consists of the variation of (t) 
with strain, for different orientations of the uniaxial deformation. In 
the last panel on the bottom right we present t/c/to and JJcj {t) for 
the two representative directions. 



FIG. 11: (Color online) Spatial variation of the magnetization as 
function of strain, (a) s = 0; (b) and (c) corresponds to s = 
0.14,0.22, with stress along the armchair edge; (d), (e), (f) corre- 
sponds to £ = 0.10, 0.14, 0.22, with stress along the zig-zag edge. 
In both cases there is an increase of the magnetization along the arm- 
chair edge as s increases. The upright triangles represent positive 
magnetization and the down ones represent negative values. 



the value of the (uniform) hopping t, a change in t produces 
a change in the absolute value of Uc. If we had an uniform 
magnetic ground state (as would be the case for bulk graphene 
without zig-zag edges), the effect of external stress could be 
captured through the average hopping (t), which diminishes 
as s increases, causing a reduction of the critical value of the 
Hubbard interaction, Uc{e) ~ (^{t{e)) [a should be around 
2.23, as per eq. ([17])]. Since we measure energies in units of 
the bare t, the above can be written as 



Uc{s) 



~ a 



(28) 



and thus Uc{e)/t would be expected to follow the variation of 
{t) with strain. 

To verify to which extent such effects contribute to the re- 
sults shown in Fig. [8] we have calculated the critical Hubbard 
interaction expected for a uniform graphene system, as a func- 
tion of magnitude and direction of strain. Within the Hartree- 
Fock framework Uc is given b^^ 



1 _ 1 ^ 1 



(29) 



where N is the total number of carbon atoms, and E{k) the 
non-interacting electron dispersion. In the presence of strain 
we will have a generalized dispersion given by 



E{k) = ± 1^2+^3 6" 



ik.ai 



ti e' 



ik.a2 I 



(30) 



reflecting that the nearest neighbor hoppings, ti, can all be 
different in general. Using the parameterization introduced 
in reference |47] we have extracted Uc as a function of strain 
magnitude, £, and orientation with respect to the honeycomb 
lattice, 0(0 = for ZZ, and (9 = 7r/2 for AC). 



Figures 10 a,b) show how the three hopping integrals ti, t2 
and ts vary under uniaxial strain along the ZZ and AC di- 
rections, respectively. Also included in those panels, are the 
respective average values of the hopping, defined as the arith- 
metic mean of the three nearest neighbor hopping integrals. 
It can be seen that the average hopping, (t), is essentially the 
same for both ZZ and AC. This is shown more clearly in 
Fig. [Toj^c), where we plot {t) as a function of strain, and for 
different strain directions: there is no sensible modification of 
(t) as the angle defining the tension direction is changed. 
Notwithstanding, the tendency is for (t) to decrease, as we 
naturally expect. The critical values of Uc under strain are 
shown in Fig. 10 ^d), where we plot both Uc/to (that reflects 
the absolute variation in the critical coupling), and Uc/{t) 
(which reflects the statement in eq. ([28])). On the one hand, 
the fact that Uc/ (t) is roughly constant up to deformations of 
20% tallies with the assumption in eq. ([28]) using a constant 
parameter a. However, even though the decrease in Uc/to 



with strain shown in Fig. 10 ^d) is qualitatively in agreement 
with the discussion above regarding the behavior of Uc for 
the dots in Fig. [9| the curves in Fig. [TOj^d) do not decrease as 
rapidly. Hence, the above argument that the critical U should 
follow the variation of (t) ([28]), is not quantitatively accurate. 

The reason for this lies in the very nature of the magnetic 
ground states of the quantum dots, which are not uniform. 
Consequently, the above argument fails in quantitative accu- 
racy, because it assumes uniformity. Similarly to what hap- 
pens in nanoribbons, the magnetization in the dots studied 
here has a strong space dependence, being highly enhanced 
near certain edges. This is a direct consequence of the charac- 
ter of the electron states around the Fermi energy, which tend 
to be localized near the boundaries. Moreover, as the dots are 
deformed by the applied strain, this space distribution is af- 



fected as well. In Fig. 11 we show the spatial evolution of 
the magnetization at the edges of the dots of type SQR2 as 
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e increases. We have chosen representative values of e such 
that the effect is clearly evident. Stress along either the arm- 
chair edge (panels (b) and (c) in Fig. [TT]), or the zig-zag edge 
(panels (d), (e) and (f) in Fig. [TT]) shows the same trend: a 
tendency for magnetization transfer from the zig-zag to the 
armchair edge of the dot. The effect is more pronounced for 
tension along the ZZ as shown in Figs.[TTJd-f). In all cases, 
the value of the magnetization increases, as can be seen in the 
lower panels of Fig. |9] From the picture described in Fig. [7] 
the behavior of the magnetization in the ZZ case can be un- 
derstood as follows: large stress along zig-zag edges tends to 
produce quasi-dimmers, weakly connected to each other; the 
dimers at the armchair edges are only coupled to the bulk of 
the dot by two weak bonds, and this favors the stabilization of 
a magnetic state. 



V. DISCUSSION 

In this paper we have studied magnetism in graphene quan- 
tum dots of particular geometries, with and without zig-zag 
edges. Similarly to what has been proposed in the context of 
graphene nanoribbons, the magnetism displayed by these sys- 
tems may be used in spin filters. We have not studied the effect 
of leads on the ferromagnetic order, but the broadening of the 
quantum dot levels due to coupling to the leads (a system with 
a continuous spectrum), will have an impact on the magnetic 
properties of the dot. In the case of a ferromagnetic nanopar- 
ticle, a theoretical description has already been developed^^, 
but no such model exists for graphene to our best knowledge. 
This study will be the subject of a forthcoming publication. 
Another aspect we have not considered in our study is the 
effect of the substrate on the magnetic properties of the dot. 
Since graphene dots are meant to be used as nanoelectronic 
devices, they will always interact with some substrate, which 
can reduce thermal fluctuations, and favor magnetism at finite 
temperatures. 

It is worth stressing again the main reason for magnetism 
at the edges of some graphene systems. Thinking about 
graphene ribbons or dots as a bulk system, that is by looking at 
the total density of states, leads immediately to the objection 
that magnetism should not be present for any finite U value. 
However, by looking at our results for magnetism in these sys- 
tems it is clear that this is a property of the edges. Therefore 
the relevant quantity is not the bulk density of states, but rather 
the local density of states, and this latter quantity, near the 
edges, does become very large at the Fermi energy {E = 0), 
thus leading to very small critical U values (eventually indis- 
tinguishable from zero). This shows that the total density of 
states is not a relevant quantity for this problem. 

Ab-initio calculations have shown that very small benzene- 
based systems, such as bisanthrene, have ferromagnetic 
ground states, and therefore we expect that magnetism will 
also be present in small quantum dots of graphene, as hinted 
by our Hartree-Fock results. We have also seen how stress 
might affect the magnetic ground states. When applied along 
the zig-zag edges, stress seems to promote a spatial rearrange- 
ment the magnetization distribution throughout the dot. This, 



combined with the transport response of these systems, may 
allow mechanical control over spin-polarized currents flow- 
ing inside the dot. On the other hand, the fact that stress along 
zig-zag edges leads to the formation of dimers weakly coupled 
between them, suggests that the system may prefer to form a 
sort of spin liquid. Investigations of whether the ground state 
will be truly magnetic or a spin liquid are required. 

In our calculations including strain, we have resorted to lin- 
ear elasticity to describe the lattice deformations, as reported 
in reference 47 We now wish to justify its validity, since it 
is an important point an deserves a careful explanation. We 
do not expect linear theory to remain valid for deformations 
as high as the 20% used in some of our calculations above. 
However, since we are interested solely in the influence of 
strain in the electronic structure, and not on the detailed elas- 
tic response, the relevant detail is not the validity of the linear 
theory itself, but rather how a certain amount of strain changes 
the hoping values. In our calculations we combined linear the- 
ory with the widespread parameterization of the dependence 
of the VppTT on the carbon-carbon distance 

Vpp^{l)=te-^-^'^'-'\ (31) 

where / is the length of the stretched (or compressed) carbon- 
carbon distance, in units of that undeformed distance. From 
the above equation what matters really is the quantity — 1), 
which is determined by the amount of strain, and its prefactor 
in the argument of the exponential. Since there might be le- 
gitimate doubts with respect to the use of linear elasticity for 
determining / in the above equation, this issue was addressed 
in another publication^. There, first principles calculations 
were used to study the effect of strain on the nearest-neighbor 
hopping values. By their nature, first principles calculations 
make no use of any elastic approximation, being valid for 
arbitrary deformations, in principle. It was found that the 
combined use of linear theory and of the above formula for 
VppTT gives accurate results for the hoping values. Moreover, 
as discussed in reference l47l subsequent ab-initio calculations 
have shown good quantitative agreement with the use of linear 
elasticity combined with the parameterization ( [3T] ) (for exam- 
ple the merging of the two Dirac points for tension along the 
zigzag direction is predicted to occur at the same values of 
strain, around 25%, both ab-initio and within tight-binding 
with linear elasticity). Therefore our approach to the calcula- 
tion of the hopping variation upon strain, and its consequences 
for the TT-bandstructure is accurately captured by the linear 
theory. 

The impact of edge roughness on the magnetic properties 
of the dots also deserves some considerations. In graphene 
nanostructures prepared with current fabrication techniques, 
the edges are invariably rough and disordered. This can hin- 
der the stability of long range magnetic order. However much 
is still unknown with respect to the nature of edge reconstruc- 
tion in real graphene nanostructures. It is expected that the 
chemical bonds at the edges be under (surface) tension, which 
leads to edge reconstruction as a means of relieving elastic 
energy. This reconstruction could lead to self-organization 
of the edge, reducing the roughness and allowing for long 
range magnetic order. Moreover, the recent advances in tai- 
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loring, at the atomic level, sub-nanostructures by transmis- 
sion electron microscope^ might allow unprecedented con- 
trol over the edge profiles. Having such control will certainly 
add greater tangibility to the prospect of tailoring magnetic 
states in graphene nanostructures. 

Finally, we point out that the experimental observation 
of ferromagnetism arising from grain boundaries with zig- 
zag edges in graphite strongly suggests that dimensionality 
is not a paramount issue in carbon-based systems, as long 
as some mechanism for quenching the thermal and quantum 
fluctuations is present. Besides the anisotropy in the inter- 
action between the local magnetic moments^^, the fact that 
graphene is generally deposited onto a substrate could pro- 
vide an extra quenching mechanism, in the same way that 



it suppresses the fluctuation-induced crumpling of the two- 
dimensional graphene membrane. The mechanical stability 
of graphene and the robustness of its magnetic phases can be 
seen as a vivid example of fluctuation quenching. 
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